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Abstract 

In this paper, an optimal batch estimator and filter based on the Minimum 
Model Error (MME) approach is developed for three-axis stabilized spacecraft. 
Three different MME algorithms are developed. The first algorithm estimates the 
attitude of a spacecraft using rate measurements. The second algorithm estimates 
the attitude without using rate measurements. The absence of rate data may be a 
result of intentional design or from unexpected failure of existing gyros. The third 
algorithm determines input-torque modeling error trajectories. All of the 
algorithms developed in this paper use attitude sensors (e.g., three-axis 
magnetometers, sun sensors, star trackers, etc). Results using these new 
algorithms indicate that an MME-based approach accurately estimates the attitude, 
rate, and input torque trajectories of an actual spacecraft. 


Introduction 

The attitude of a spacecraft can be determined by either deterministic methods or by utilizing 
algorithms which combine dynamic models with sensor data. Three-axis deterministic methods, such as 
TRIAD [1], QUEST [2], and FOAM [3], require at least two simultaneous vector measurements to 
determine the attitude (direction-cosine) matrix. An advantage of both the QUEST and FOAM 
algorithms is that the attitude of a spacecraft can be estimated using more than two measurements. This 
is accomplished by minimizing a quadratic loss function first posed by Wahba [4]. However, all 
deterministic methods fail when only one vector measurement is available, (e.g., magnetometer data 
only). Estimation algorithms utilize dynamic models, and subsequently can (in theory) estimate the 
attitude of a spacecraft using measurements of a single reference vector. Although all spacecraft in use 
today have at least two on-board attitude sensors, estimation techniques can be used to determine the 
attitude during anomalous periods, such as solar eclipse and/or sensor co-alignment 

The most commonly used technique for attitude estimation is the Kalman filter [5]. The Kalman filter 
utilizes state-space representations to both estimate plant dynamics and also filter noisy data. Errors in 
the dynamical model and measurement process are assumed to be modeled by a zero-mean Gaussian 
process with known covariance. The optimality criterion in the Kalman filter minimizes the trace error 
covariance between estimated responses and model responses. In theory, the Kalman filter does not 
require actual measurements to satisfy this optimality criterion; however, in actual practice measurements 
are often used to properly “tune” the filter estimates. 

Smoothing algorithms further refine state estimates by utilizing both a “forward filter” and a 
“backward filter” (see e.g., Gelb [6]). An advantage of smoothing algorithms is that the error covariance 
is always less or equal to either the forward or backward filter alone. A disadvantage of smoothing 
algorithms is that they cannot be implemented in sequential (real-time) estimation. 

In order for the Kalman filter to be truly optimal, both the measurement error process and model 
error process must be random Gaussian processes with known covariance. In most circumstances. 
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properties of the measurement error process are known a priori by utilizing statistical inferences applied 
to sensor measurements. However, model error statistics are not usually weU known. In actual practice 
the determination of the model error covariance in the Kalman filter is usually obtained by an ad hoc 
and/or heuristic approach, which can result in suboptimal filter designs (e.g., determining random gyro 
drift rate). Also, in many instances, such as nonlinear model errors or non-stationary processes, the 
assumption of a stationary Gaussian process can lead to severely degraded state estimates. 

For spacecraft attitude estimation, the Kalman filter is most applicable to spacecraft equipped with 
three-axis gyros as well as attitude sensors [7]. However, rate gyros are generally expensive and are 
often prone to degradation or failure. Therefore, in recent years rate gyros have been omitted (e.g., in 
Small Explorer (SMEX) spacecraft, such as Solar Anomalous Magnetospheric Particle Explorer 
(S AMPEX) spacecraft). To circumvent the problem of rate gyro omission or failure, analytical models of 
rate motion can be used. This approach has been successfully used in a Real-Time Sequential Filter 
(RTSF) algorithm which propagates state estimates and error covariances using dynamic models [8]. The 
estimation of dynamic rates by the RTSF is accomplished from angular momentum model propagation, 
and then correcting for these rates by using a “gyro bias” component in the filter design. A clear 
advantage of using dynamic models is shown for the case of sun-magnetic field near co-alignment. For 
this case, deterministic algorithms, such as TRIAD and QUEST, show anomalous behaviors with extreme 
deviations in determined attitudes. Since the RTSF propagates an analytical model of motion, attitude 
estimates are improved even when data from only one attitude sensor is available. However, the RTSF is 
essentially a Kalman filter in which the “gyro bias” model (and subsequently the angular momentum 
model correction) is assumed to be a Gaussian process with known covariance. Also, fairly accurate 
models of angular momentum are required in order to obtain accurate estimates. Subsequently, the 
design process for choosing the model error covariance becomes difficult 

In this paper, an optimal attitude estimation algorithm is developed which is capable of robust and 
accurate state estimation for spacecraft lacking accurate or any rate measurements and/or accurate 
dynamic models. This algorithm is based on the Minimum Model Error (MME) [9] batch-estimation 
approach. The advantages of the MME estimator over conventional Kalman strategies include: (i) no a 
priori statistics on the form of the model error are required, (ii) the actual model error is determined as 
part of the solution, and (iii) state estimates are free of jump discontinuities, which greatly smoothes out 
high measurement noise. The MME estimation approach has been successfully applied to numerous 
poorly-modeled dynamic systems which exhibit highly nonlinear behaviors (see, e.g. [10-11]). Previous 
MME studies used TRIAD-determined quaternions as measurements [12]. The formulations developed 
in this paper expand upon this method to include attitude sensors, such as three-axis magnetometers 
(TAM), fine sun sensors (FSS), star trackers, etc. 

The organization of this paper proceeds as follows. First, a summary of the spacecraft attitude 
kinematics and sensor models is shown. Then, a brief review of the MME estimator for nonlinear 
systems is shown. Next, various MME-based algorithms are developed for the purpose of attitude 
estimation, which include: a simple linear algorithm which is used to smooth noisy rate measurements, an 
attitude estimator using rate measurement information, an attitude estimator without the utilization of any 
rate meas urements, and an input torque estimator. Finally, these MME designs are used to estimate the 
attitude, rate, and input torque trajectories of the SAMPEX spacecraft in order to demonstrate the 
usefulness of these algorithms. 


Attitude Kinematics and Dynamics 

In this section, a brief review of the kinematic and dynamic equations of motion for a three-axis 
stabilized spacecraft is shown. The attitude is assumed to be represented by a quaternion, defined as 
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where n is a unit vector corresponding to the axis of rotation and 0 is the angle of rotation. The 
quaternion kinematic equations of motion are derived by using the spacecraft’s angular velocity (co), 
given by ~ 

? = = (3) 


where ti(co) and e(^) are defined as 
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The 3x3 dimensional matrices [ax] and [^ ]3 X ] are referred to as cross product matrices since 
axb = [ax]b, with 
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Since a three degree-of-freedom attitude system is represented by a four-dimensional vector, the 
quaternions cannot be independent. This condition leads to the following normalization constraint 


Q 7 ! = ^13^13 + *4=1 

( 6 ) 

Also, the matrix obeys the following helpful relations 
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The dynamic equations of motion, also known as Euler s equations, for a rotating spacecraft are 
given by ([13]) 

L = N-co xL = /^cp (8) 

where L is the total angular momentum, N_ is the total external torque (which includes, e.g., control 
torques, aerodynamic drag torques, solar pressure torques, etc.), and l b is the inertia matrix of the 
spacecraft. If reaction or momentum wheels are used on the spacecraft, the total angular momentum is 
given by 

Ll = + 

where h is the total angular momentum due to the wheels. Thus, Equation (8) can be re-wntten as 

L = K-[ib l (L-h)\ x L ( 10 > 

The measurement model is assumed to be of the form given by 

B B = A(q)B, (11) 

where is a 3x1 dimensional vector of some reference object (e.g., a vector to the sun or to a star, or 
the Earth’s magnetic field vector) in a reference coordinate system, B b is a 3x1 dimensional vector 
defining the components of the corresponding reference vector measured in the spacecraft body frame, 
and a(<?) is given by 

= (<& ^13— 13 ) ^3x3 + 2 ^i 3 ^[ 3 ~ 2<?4 fel3 X ] (12) 

which is the 3x3 dimensional (orthogonal) attitude matrix. 


Minimum Model Error Estimation 

In this section, a brief review of the Minimum Model Error (MME) estimation algorithm is shown. 
The essential feature of this batch estimator is that actual model error trajectories are determined during 
the estimation process, unlike most filter/smoother algorithms which assume that the model error is a 
stochastic process with known properties. The MME algorithm determines the correction added to the 
assumed model which yields an accurate representation of the system’s behavior. This is accomplished by 
solving an op timali ty condition using an output residual constraint. Therefore, accurate state estimates 
can be determined without the use of precise system representations in the assumed model. 

The MME algorithm assumes that the state estimates are given by a preliminary model and a to-be- 
determined model error vector, given by 

1(0 = /[l(0> «(0> d{t\ r] ( 1 3a) 

K0 = £[*(')>'] (13b) 


where / is an nxl model vector, x(t) is an nxl state estimate vector, u(t) is a pxl vector of known 
inputs, and d(t) is an nxl model error vector, £ is a x l measurement (sensitivity) vector, and y(t) is a 
q x\ estimated output vector. State-observable discrete measurements are assumed for Equation (13b) in 
the following form 
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(14) 


l{tk) = 8 k [x(tk)>tk]+Vk 

where y(q) is a qxl measurement vector at time t k , and v k is a qxl measurement noise vector which is 

assumed to be a zero-mean, Gaussian distributed process with known covariance. 

In the MME algorithm, die optimal state estimates are determined on the basis that the measurement- 
minus-estimate error covariance matrix must match the measurement-minus-truth error covariance 
matrix. This condition is referred to as the “covariance constraint,” shown as 

= R k (15) 

where R k is the element-by-element (known) measurement error covariance. However, problems may 
arise using Equation (15) which are attributed to “small sample” statistics [14]. Therefore, in the typical 
case where the measurement error process is stationary, the average covariance can be used, given by 

1 m 
*=i 


where m is the total number of measurements. 

Next, the following cost function is minimized with respect to d(x) 

1 m T l f 

R l {£('*)"* \d T {x)Wd{x)dx (17) 

k=\ 1 : 

l 0 

where W is an nxn positive-definite weighting matrix. The necessary conditions for the minimization of 
Equation (17) lead to the following two-point-boundary-value-problem (TPBVP) [9] 

1(0 = /U(0> «(0. d(t), t ] ( 1 8a) 
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where X(t) is an nxi co-state vector which is updated at each measurement point using Equation (18d). 
The boundary conditions are selected such that either X(»o) = 0 or i(r 0 ) is specified at the initial time and 

either \(t } ) = 0 or x(t f ) is specified at the final time. 

The solution of the TPB VP for a given weighting matrix yields a state estimate time trajectory which 
can be used to determine a measurement residual covariance matrix. The covariance constraint is 
satisfied when the proper balance between model error and measurement residual has been achieved. If 
the measurement residual covariance is higher than the known measurement error covariance (R), then IV 
should be decreased to less penalize the model error. Conversely, if the residual covariance is lower than 
the known covariance, then W should be increased so that less unmodeled dynamics are added to the 
assumed system model. The optimal weighting matrix is therefore obtained when the covariance 
constraint in Equation (16) is satisfied. 


Gvro Noise Smoother 

Gyros tend to be noisy and have an inherent drift. Also, gyros are usually sampled at a higher 
frequency than attitude sensors. In order to first filter the noise, a simple MME-based smoothing 
algorithm can first be applied. This algorithm minimizes 

m tf 
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where & g (t) is the estimated gyro output, and d g (t) is the model error correction. Minimizing Equation 
(19) leads to the following TPBVP 
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The solution of Equation (21) can be determined by using a steady-state Riccati transformation (see [15] 
for details). This transformation leads to the following 
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where the subscript (i) represents each gyro measurement set, and A/ is the sampling interval. For a 
given weighting and measurement covariance, the first step is to determine the steady-state Riccati 
solution in using Equation (22a). Then, the inhomogeneous Riccati trajectory is solved backwards in 
time using Equation (22b), with discrete jumps at each measurement point given by Equation (22c). 
Finally, the smoothed gyro estimates are determined using Equation (22d). An advantage of this 
algorithm is not only the inherent smoothing properties, but also that the gyro estimates are totally 
continuous. Therefore, the generally discrete gyro measurements can be replaced with the continuous 
gyro estimates given by Equation (22d). 


Attitude Estimation Using Rate Measurements 

The MME attitude angular velocity estimation formulation using rate measurements minimizes the 
following cost function 


1 m T */ 

J = 2 -4?)S/} +±\d T (i)Wd(T)dx 

k= 1 '* tk , 


(23) 


subject to 

|(0 = \ g (0 + m] £(0. = £ 0 


(24) 


where Gi g (t) is the rate measurement vector, q{t) is the estimated quaternion, and B b and B, are the 

spacecraft body measurement and corresponding inertial field vector, respectively. The model error (d) 
is a correction to the rate measurements, which forces the model responses to satisfy the covariance 
constraint in Equation (16). 

The TPBVP given by Equation (18) can be written as 
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The sensitivity matrix H in Equation (25d) can be derived as 
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(27a) 





-^ 3x3 + [^ 13 x ] 



(27b) 


The extension to using multiple attitude sensors is accomplished by using a partitioned residual output 
and sensitivity matrix, given by 



- 41 ) 8 ,,} 


(28) 


The co-state update in Equation (25d) shows a nonlinear relationship with respect to the quaternion 
estimate. However, this nonlinearity can be reduced to be a linear function if the quaternions obey 
normalization and each attitude sensor is assumed isotropic. This can be shown by deriving the co-state 
update using 


_i__a_ 

2 r dq 




(29) 


where the measurement covariance is now assumed to be isotropic for each sensor (i.e., the measurement 
errors in each one of the axes are assumed equal). Therefore, /? = r/ 3x3 , which is a valid assumption for 
almost all attitude sensors. In order to determine the partial derivative in Equation (29), the following 
identities and definitions are used 
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Equation (29) can now be re-written as 
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The partial derivative in Equation (31) is given by 
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Hence, if the quaternions obey normalization the following identity is true 

{n(I s )E(fl ; )|+(sf b,)|} = m{i B -*(«)*,} 

Therefore, if the sensor measurements are isotropic, the co-state update in Equation (25d) is linear with 
respect to the quaternion estimate. 

"Die TPBVP shown in Equations (25a)- (25d) can be solved by using gradient techniques. The basic 
gradient procedure is to first guess for the model error trajectory (d). Then, integrate the quaternion 
states forward using Equation (25a) and co-states backward using Equation (25c) accounting for discrete 
jumps in Equation (25d). The next search direction is given by Equation (25b). This procedure is 
continued until convergence is achieved. 


Attitude Estimation without Rate Measurements 

In this section, the MME estimator is derived for spacecraft which lack any rate information. The 
formulation is based upon using Euler’s equation for modeling the angular momentum. The MME 
problem for this case minimizes the following cost function 
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where 
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(36) 


where h_ is the measured angular momentum due to the wheels. Minimizing Equation (34) leads to the 
following TPBVP 
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with discrete jumps in the co-states given by 
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(38) 


M'* )■ - M'* ) + H \ *''{£« - M&iX'- M'/)=Q 

The TPBVP given by Equations (37) and (38) can be solved by using a simple gradient-based search 
technique. 


Input Torque Estimation 

In this section, the MME estimator is used to estimate model error torques using angular rate 
trajectories. These angular rate trajectories are assumed to be known (e.g., from finite differenced 
attitude es tima tes, or from angular rate estimates from an MME design or other estimator). First, a 
measured angular momentum vector is determined by 

L = IbQ+lL (39) 


In general, the angular momentum measurements in Equation (39) will be noisy due to the angular 
velocity measurements of the wheel speed. However, this noise is inherently smoothed by the MME 
estimator. The MME problem for determining the errors in the torque input of Euler’s equation 
minimiz es the following cost function 


J = -^{L-L} T R- 1 {L-L\ +^jd T (x)Wd(x)dx 

k = 1 tk * /n 


(40) 


subject to 
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Mi nimizin g Equation (40) leads to the following TPBVP 
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The solution to the TPBVP in Equation (42) can be determined by using a Riccati transformation [15]. 
Applying this technique leads to the following equations 
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hi) = {-[<0(0 x]- W~' />(!) }£(()- »'“' 4(<)+M<), U‘o)=Lo (43e) 

Therefore, the first step is to solve for the Riccati and inhomogeneous trajectories backwards in time 
using Equations (43a) and (43c), accounting for discrete jumps by Equations (43b) and (43d). Then, the 
angular momentum estimates are determined by integrating Equation (43e) forwards in time. 


Attitude Estimation of an Actual Spacecraft 

In this section, the MME estimation algorithms previously developed are used to estimate the 
attitude, rate, and input torque trajectories of the SAMPEX spacecraft using vector measurement 
observations. The SAMPEX general mission is to study energetic particles and various types of rays. 
The spacecraft is three-axis stabilized in a 550 by 675 km elliptical orbit with an 82° inclination. The 
attitude control hardware consists of a magnetic torquer assembly (MTA) and a reaction wheel assembly 
(RWA). The attitude determination hardware consists of five coarse Sun sensors (CSS) (primarily for 
Sun-acquisition), one fine Sun sensor (FSS), and a three-axis magnetometer (TAM). Also, no rate 
gyroscopic instruments are present on the spacecraft. 

The onboard computer routine to determine attitude is based upon the TRIAD [1] deterministic 
method. The spacecraft is controlled by the MTA to maintain the fixed solar arrays perpendicular to the 
sun-line. The RWA is used to point the instrument boresight axis as required by the scientific mission. 
During eclipse no sun measurements are available from the FSS. Attitude control is maintained by using 
a constant sun-line vector as a “pseudo-measurement,” so that both the MTA and RWA are still utilized. 
During vector co-alignment, the spacecraft is placed in a “coast” mode in which the MTA is not used (see 
[16] for more details). The required nominal attitude determination accuracy is ±2°. During anomalous 
conditions (eclipse and/or measurement vector co-alignment) the attitude cannot be determined by 
deterministic methods, such as TRIAD. The MME algorithms presented in this paper can determine the 
attitude using TAM measurements only, so that attitude accuracy may be checked for any deviations from 
nominal performance. 

The inertial field trajectories are obtained by using a 8th order spherical harmonic model of the 
Earth’s magnetic field with International Geomagnetic Reference Field (IGRF) coefficients. 
Magnetometer measurements by the TAM are known to be extremely accurate (within 0.3 mG). 
However, experience has shown that errors in the magnetic field model have a standard deviation of 
about 3 mG [17]. Therefore, 9 mG 2 is chosen for the diagonal elements of the measurement covariance 
matrix. 

The first test case involves using both TAM and FSS measurements. A plot of the finite differenced 
angular rates using TRIAD determined attitudes is shown in Figure 1. These rates are extremely noisy, 
which is due to the large digitization noise associated with the FSS measurements. The TRIAD 
determined rates are next used in the first MME formulation (with rate information), along with the TAM 
and FSS measurements. A plot of the MME estimated rates is shown in Figure 2. Clearly, these rates are 
smoother than the TRIAD determined rates. Next, the MME input torque estimator is applied using 
these estimated rates. A plot of the MME determined input torques is shown in Figure 3. These torques 
correspond to a correction to the dynamic model, so that the model responses match the vector 
measurement observations. 

The second case involves using TAM measurements only to estimate the attitude and angular rates. 
The MME attitude estimator without rate measurements is used for this case. A plot of the estimated 
angular rate trajectories is shown in Figure 4. These angular rate estimates clearly show a rotation about 
the spacecraft’s y-axis, which is the desired motion. A plot of the error between the estimated MME 
attitudes and the attitudes determined by TRIAD is shown in Figure 5. A slight hangoff is seen in the 
pitch axis. This may be due to nonlinear effects in the magnetic field model (this hangoff is also seen in 
Kalman filter approaches for other spacecraft such as UARS). However, the MME algorithm is able to 
determine attitudes to within 0.3° using TAM data only. This can be useful in determining the attitude 
when deterministic methods fail. 
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Conclusions 

In this paper, several MME algorithms were presented for use in attitude estimation using vector 
measurement observations. The first algorithm used angular rate measurements to determine attitude 
estimate trajectories. The second algorithm estimated the attitude trajectories without any rate 
measurement information. The third algorithm determined the required torque input trajectories so that 
the model responses match the vector observations. An advantage of all of these algorithms is that 
quaternion normalization was maintained, since linearization of the dynamic model was not needed. The 
MME-based algorithms were then applied to an actual spacecraft. Results indicated that an MME-based 
approach provides a robust algorithm which can be used to determine the attitude, rate, and modeling 
error torque trajectories of a spacecraft from vector measurements. 
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Figure 1 Plot of TRIAD Determined Attitude Rates 
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Figure 4 Plot of MME Estimated Rates Using TAM Data Only 
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Figure 5 Plot of Attitude Errors Between TRIAD and MME 
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